# for tidy style coding and plottinglibrary(tidyverse) library(vroom) # to read lots of csv files at once# more table optionslibrary(pander) # for tableslibrary(kableExtra) # for scrolling tableslibrary(data.table) # for efficient handling of large dataframes# making ggplot more powerfullibrary(MetBrewer) # for colour palettes based upon artwork housed at the METlibrary(MoMAColors) # for colour palettes based upon artwork housed at MoMAlibrary(wesanderson) # for colour palettes based on wes anderson movieslibrary(rcartocolor) # for nice sequential colour schemeslibrary(PNWColors) # for colour palettes library(tidybayes) # for plotting distributionslibrary(stickylabeller) # labelling facets with strings in ggplotlibrary(geomtextpath) # for curved plot annotationslibrary(ggtext) # for markdown syntax in plot labelslibrary(patchwork) # for patching plots togetherlibrary(ggnewscale) # to reset scales in plots, allowing multiple fill arguments in ggplot# for computation speed checkslibrary(profvis) # breakdown of complex functionslibrary(bench) # individual functions
Plotting results from the inclusive fitness model
We find the inbreeding depression tolerance thresholds with strict polygyny are
where \(r\) is the genealogical relatedness coefficient, \(c_{\mathrm{m},i}\) and \(c_{\mathrm{f},i}\) are the juvenile class reproductive values of males and females, respectively, and \(F_\mathrm{f}\) and \(F_\mathrm{m}\) are the female and male inbreeding coefficients. See main text for more details.
In the bottom panel of Figure 1, we show the discrepancy between the sexes, for what level of inbreeding depression can be tolerated before inbreeding acceptance alleles are not favoured by selection.
Code
sexual_conflict_plot <- figure_1_thresholds %>%filter(r ==0.5) %>%mutate(Case =case_when(str_detect(Case, "W") ~"W",str_detect(Case, "Y") ~"Y",.default = Case),Sex =case_when(str_detect(Case, "Cyto") ~"Both",str_detect(Sex,"Female") ~"Female",str_detect(Sex, "Male") ~"Male"),D =case_when(str_detect(Case, "p") ~ D -0.000001,.default = D)) %>%ggplot(aes(x =1-D, y =reorder(Case, -D), group = Case)) +geom_line(size =1) +geom_point(aes(fill = Sex), shape =21, size =9, stroke =1) +scale_fill_manual(values =c("#d0e2af", "pink", "aliceblue")) +annotate("text", x =0, y =1, label ="0") +annotate("text", x =0, y =2, label ="0") +annotate("text", x =0.2, y =3, label ="1/5") +annotate("text", x =0.5, y =3, label ="1/2") +annotate("text", x =0.33, y =4, label ="1/3") +annotate("text", x =0.67, y =4, label ="2/3") +annotate("text", x =0.5, y =5, label ="1/2") +annotate("text", x =0.8, y =5, label ="4/5") +annotate("text", x =1, y =6, label ="1") +annotate("text", x =1, y =7, label ="1") +labs(x =~delta~'(inbreeding depression)',y ="Genomic location",title ="C . Sexual conflict over full-sibling inbreeding") +scale_x_continuous(limits =c(-0.05, 1.05)) +theme_bw() +theme(axis.title =element_text(size =16),axis.text =element_text(size =12),legend.position ="bottom",legend.title =element_text(size =12),legend.text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14))sexual_conflict_plot
Figure 1. Within-organism conflict over inbreeding acceptance. In panels A and B, contours show the the severity of inbreeding depression below which inbreeding increases the inclusive fitness of genes in a particular genomic region. As the fill gets darker, fewer parts of the genome have higher fitness when inbreeding is preferred over outcrossing. C shows the range of inbreeding depression intensities that produce sexual conflict over inbreeding acceptance, when mating is between full-siblings (\(r=0.5\)). Lines indicate the conflict zone when loci affecting female and male trait values are found within the same genomic region. Values within the coloured circles are \(\delta_{\mathrm{f},i}^*\) and \(\delta_{\mathrm{f},i}^*\). Green points show where the sexes have aligned fitness interests. Results are shown for the frequent inbreeding scenario. Cyto. (m) and Cyto. (p) denote cytoplasmic genes that have strict maternal and paternal inheritance, respectively.
The simulation
Build convenience functions
A sampling function that can handle vectors of length one
Code
# so we can sample from vectors with length 1 without this being interpreted as an integersample_vec <-function(x, ...) x[sample(length(x), ...)]
A function that builds the appropriate inheritance system.
Code
make_mating_table <-function(gene_location){ make_offspring <-function(X, Y, offspring_genotype, zygote_freq, sex){data.frame(Female_genotype = X,Male_genotype = Y, offspring_genotype, zygote_freq, sex) }# Specify the possible offspring genotypes for all the potential crosses; we use these for the offspring_genotype argument in the make_offspring function# offspring genotypes offspring_genotypes_1 <-c(2,2) offspring_genotypes_2 <-c(1, 1, 2, 2) offspring_genotypes_3 <-c(1, 1) offspring_genotypes_4 <-c(0, 0, 1, 1, 2, 2) offspring_genotypes_5 <-c(0, 0, 1, 1) offspring_genotypes_6 <-c(0, 0) offspring_genotypes_7 <-c(1, 2) offspring_genotypes_8 <-c(0, 1, 1, 2) offspring_genotypes_9 <-c(0, 1) offspring_genotypes_10 <-c(1, 0) # this is diff from above bc of the order with the sexes offspring_genotypes_11 <-c(2, 1) offspring_genotypes_12 <-c(2,1,1,0)# offspring sex offspring_sex_2 <-c(0, 1) offspring_sex_4 <-c(0, 1, 0, 1) offspring_sex_6 <-c(0, 1, 0, 1, 0, 1)# even frequency of two offspring genotypes freq_2 <-rep(0.5, 2)# even frequency between four offspring types freq_4 <-rep(0.25, 4)# when there are 6 offspring genotypes freq_6 <-c(0.125, 0.125,0.25, 0.25,0.125, 0.125)if(gene_location =="A"){ books <-rbind(make_offspring(2, 2, offspring_genotypes_1, freq_2, offspring_sex_2),make_offspring(2, 1, offspring_genotypes_2, freq_4, offspring_sex_4),make_offspring(2, 0, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(1, 2, offspring_genotypes_2, freq_4, offspring_sex_4),make_offspring(1, 1, offspring_genotypes_4, freq_6, offspring_sex_6),make_offspring(1, 0, offspring_genotypes_5, freq_4, offspring_sex_4),make_offspring(0, 2, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(0, 1, offspring_genotypes_5, freq_4, offspring_sex_4),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="X"){ books <-rbind(make_offspring(2, 1, offspring_genotypes_7, freq_2, offspring_sex_2),make_offspring(2, 0, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(1, 1, offspring_genotypes_8, freq_4, offspring_sex_4),make_offspring(1, 0, offspring_genotypes_5, freq_4, offspring_sex_4),make_offspring(0, 1, offspring_genotypes_9, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="Y"){ books <-rbind(make_offspring(0, 1, offspring_genotypes_10, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="Z"){ books <-rbind(make_offspring(1, 2, offspring_genotypes_11, freq_2, offspring_sex_2),make_offspring(0, 2, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(1, 1, offspring_genotypes_12, freq_4, offspring_sex_4),make_offspring(0, 1, offspring_genotypes_5, freq_4, offspring_sex_4),make_offspring(1, 0, offspring_genotypes_10, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="W"){ books <-rbind(make_offspring(1, 0, offspring_genotypes_9, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="C"){ books <-rbind(make_offspring(1, 0, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(1, 1, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2),make_offspring(0, 1, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="P"){ books <-rbind(make_offspring(1, 0, offspring_genotypes_6, freq_2, offspring_sex_2),make_offspring(1, 1, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2),make_offspring(0, 1, offspring_genotypes_3, freq_2, offspring_sex_2) ) }return(books) }offspring_genotypes_autosome <-make_mating_table("A")offspring_genotypes_X <-make_mating_table("X")offspring_genotypes_Y <-make_mating_table("Y")offspring_genotypes_Z <-make_mating_table("Z")offspring_genotypes_W <-make_mating_table("W")offspring_genotypes_C <-make_mating_table("C")offspring_genotypes_P <-make_mating_table("P")
A function that takes two parental genotypes and produces offspring
Code
sample_mating_table <-function(inheritance_scheme, f, mother){# cut to possible genotypes possibilities <- inheritance_scheme[inheritance_scheme$Female_genotype == mother[4] & inheritance_scheme$Male_genotype == mother[9], c(3,5)]# get prob of producing each genotype probs <- inheritance_scheme[inheritance_scheme$Female_genotype == mother[4] & inheritance_scheme$Male_genotype == mother[9], 4]# sample possibilities[sample(size = f,x =nrow(possibilities), prob = probs,replace =TRUE), ]}
Load the parameter space
Code
resolution <-25starting_pop_size_autosomes <-2000# both sexes harbour two copies of each autosomal chromosome = 1000 autosomal haplotypesparameters <-expand_grid(chromosome =c("A", "X", "Y", "Z", "W", "C", "P"),v =c(8, 80),D =seq(0, -0.99, length = resolution), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.01, 1, length = resolution) ) %>%full_join(tibble(chromosome =c("A", "A", "A", "A", "A", "A","X", "X", "X", "X","Y","Z", "Z", "Z", "Z","W","C", "C","P", "P"),sex_expressed =c(0, 0, 0, 1, 1, 1,0, 1, 1, 1,0,0, 0, 0, 1,1,0, 1,0, 1),dominance =c(0, 0.5, 1, 0, 0.5, 1,1, 0, 0.5, 1, 1, 0, 0.5, 1, 1,1,1, 1,1, 1)) %>%mutate(Starting_pop_size =case_when(chromosome =="A"~ starting_pop_size_autosomes, chromosome =="X"| chromosome =="Z"~ starting_pop_size_autosomes /0.75, chromosome =="Y"| chromosome =="W"~ starting_pop_size_autosomes*4, chromosome =="C"| chromosome =="P"~ starting_pop_size_autosomes*2)),relationship ="many-to-many", by ="chromosome") %>%mutate(baseline_mean_lifespan =1,v = v / (Starting_pop_size /2),f =5, refractory_period =-log(refractory_period_prop_cohort_alive),mutation_time =5, # this is when the mutation can be introduced fromtime_end =1000, # with avg lifespan = 1, this is ~ roughly 1000 gensparameter_space_ID =row_number(),mutation_events =5)parameters_autosome <- parameters %>%filter(chromosome =="A") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_X <- parameters %>%filter(chromosome =="X") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_Y <- parameters %>%filter(chromosome =="Y") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_Z <- parameters %>%filter(chromosome =="Z") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_W <- parameters %>%filter(chromosome =="W") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_C <- parameters %>%filter(chromosome =="C") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_P <- parameters %>%filter(chromosome =="P") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsif(!file.exists("parameters/parameters_autosome.txt")){ parameters_autosome %>%write.table("parameters/parameters_autosome.txt") parameters_X %>%write.table("parameters/parameters_X.txt") parameters_Y %>%write.table("parameters/parameters_Y.txt") parameters_Z %>%write.table("parameters/parameters_Z.txt") parameters_W %>%write.table("parameters/parameters_W.txt") parameters_C %>%write.table("parameters/parameters_C.txt") parameters_P %>%write.table("parameters/parameters_P.txt")}
We also run the simulation across a more specific parameter space, that corresponds to the zone of uncertainty produced by the dynamic effect of inbreeding on relatedness.
Code
resolution_fine <-80parameters_fine_X <-bind_rows(expand_grid(chromosome ="X",v =80,D =seq(-0.5, -0.67, length = resolution_fine), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.8, 0.999, length =15) ) %>%full_join(tibble(chromosome ="X",sex_expressed =0,dominance =1),relationship ="many-to-many", by ="chromosome"),expand_grid(chromosome ="X",v =80,D =seq(-0.2, -0.33, length = resolution_fine), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.8, 0.999, length =15) ) %>%full_join(tibble(chromosome ="X",sex_expressed =1,dominance =1),relationship ="many-to-many", by ="chromosome") ) %>%mutate(Starting_pop_size = starting_pop_size_autosomes /0.75,baseline_mean_lifespan =1,v = v / (Starting_pop_size /2),f =5, refractory_period =-log(refractory_period_prop_cohort_alive),mutation_time =5, # this is when the mutation can be introduced fromtime_end =1000, # with avg lifespan = 1, this is ~ roughly 1000 gensparameter_space_ID =row_number(),mutation_events =5) %>%slice_sample(prop =1)parameters_fine_Z <-bind_rows(expand_grid(chromosome ="Z",v =80,D =seq(-0.67, -0.8, length = resolution_fine), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.8, 0.999, length =15) ) %>%full_join(tibble(chromosome ="Z",sex_expressed =0,dominance =1),relationship ="many-to-many", by ="chromosome"),expand_grid(chromosome ="Z",v =80,D =seq(-0.33, -0.5, length = resolution_fine), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.8, 0.999, length =15) ) %>%full_join(tibble(chromosome ="Z",sex_expressed =1,dominance =1),relationship ="many-to-many", by ="chromosome") ) %>%mutate(Starting_pop_size = starting_pop_size_autosomes /0.75,baseline_mean_lifespan =1,v = v / (Starting_pop_size /2),f =5, refractory_period =-log(refractory_period_prop_cohort_alive),mutation_time =5, # this is when the mutation can be introduced fromtime_end =1000, # with avg lifespan = 1, this is ~ roughly 1000 gensparameter_space_ID =row_number(),mutation_events =5) %>%slice_sample(prop =1)parameters_fine_AX <-bind_rows(expand_grid(chromosome ="A",v =80,D =seq(-0.5, -0.67, length = resolution_fine), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.8, 0.999, length =15) ) %>%full_join(tibble(chromosome ="A",sex_expressed =0,dominance =1),relationship ="many-to-many", by ="chromosome"),expand_grid(chromosome ="A",v =80,D =seq(-0.2, -0.33, length = resolution_fine), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.8, 0.999, length =15) ) %>%full_join(tibble(chromosome ="A",sex_expressed =1,dominance =1),relationship ="many-to-many", by ="chromosome") ) %>%mutate(Starting_pop_size = starting_pop_size_autosomes,baseline_mean_lifespan =1,v = v / (Starting_pop_size /2),f =5, refractory_period =-log(refractory_period_prop_cohort_alive),mutation_time =5, # this is when the mutation can be introduced fromtime_end =1000, # with avg lifespan = 1, this is ~ roughly 1000 gensparameter_space_ID =row_number(),mutation_events =5) %>%slice_sample(prop =1)parameters_fine_AZ <-bind_rows(expand_grid(chromosome ="A",v =80,D =seq(-0.67, -0.8, length = resolution_fine), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.8, 0.999, length =15) ) %>%full_join(tibble(chromosome ="A",sex_expressed =0,dominance =1),relationship ="many-to-many", by ="chromosome"),expand_grid(chromosome ="A",v =80,D =seq(-0.33, -0.5, length = resolution_fine), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.8, 0.999, length =15) ) %>%full_join(tibble(chromosome ="A",sex_expressed =1,dominance =1),relationship ="many-to-many", by ="chromosome") ) %>%mutate(Starting_pop_size = starting_pop_size_autosomes,baseline_mean_lifespan =1,v = v / (Starting_pop_size /2),f =5, refractory_period =-log(refractory_period_prop_cohort_alive),mutation_time =5, # this is when the mutation can be introduced fromtime_end =1000, # with avg lifespan = 1, this is ~ roughly 1000 gensparameter_space_ID =row_number(),mutation_events =5) %>%slice_sample(prop =1)if(!file.exists("parameters/parameters_fine_X.txt")){ parameters_fine_X %>%write.table("parameters/parameters_fine_X.txt") parameters_fine_Z %>%write.table("parameters/parameters_fine_Z.txt") parameters_fine_AX %>%write.table("parameters/parameters_fine_AX.txt") parameters_fine_AZ %>%write.table("parameters/parameters_fine_AZ.txt")}
The simulation function
Code
continuous_time_simulation <-function(row, parameters, inheritance_scheme){print(paste("Doing row", row)) # this shows which row in the parameter space is being modelled Starting_pop_size <-round(parameters$Starting_pop_size[row], 0) f <- parameters$f[row] # fecundity constant mutation_time <- parameters$mutation_time[row] # introduce an I allele after family structure is established baseline_mean_lifespan <- parameters$baseline_mean_lifespan[row] # constant at 1 time_end <- parameters$time_end[row] # a cut-off point for each run sex_expressed <- parameters$sex_expressed[row] chromosome <- parameters$chromosome[row] v <- parameters$v[row] refractory_period <- parameters$refractory_period[row] D <- parameters$D[row] dominance <- parameters$dominance[row] parameter_space_ID <- parameters$parameter_space_ID[row] mutation_events <- parameters$mutation_events[row]# Set the number of breeding sites breeding_sites <-round(0.2*Starting_pop_size, 0)# what inheritance system does this run follow offspring_genotypes <- inheritance_scheme# Set the maximum number of I alleles that can be found in each sexif(chromosome =="A"){ female_max_I <-2 male_max_I <-2 }if(chromosome =="X"){ female_max_I <-2 male_max_I <-1 }if(chromosome =="Z"){ female_max_I <-1 male_max_I <-2 }if(chromosome =="Y"){ female_max_I <-0 male_max_I <-1 }if(chromosome =="W"){ female_max_I <-1 male_max_I <-0 }if(chromosome =="C"| chromosome =="P"){ female_max_I <-1 male_max_I <-1 }# make matrix to hold results; updated as sim progresses# col1 = time, col2 = prop I, col3 = pop size, col4 = prop virgin female deaths results_matrix <-matrix(nrow = time_end*4+2, ncol =4) # record each time point# make matrix to hold population; updated as sim progresses# col1 = ID # col2 = Family ID# col3 = Sex: females = 1 and males = 0# col4 = Genotype: 0, 1 and 2 = copies of inbreeding allele# col5 = mortality rate# col6 = encountered relative: NA = NO, 1 = YES# col7 = mating state: -Inf not in pop, NA = unmated, real = out, Inf = mated female# col8 = inbred mating: NA = NO, 1 = YES (only matters for females)# col9 = mated_genotype: NA = unmated, otherwise 0,1,2 (see mating table)# col10 = breeding site: NA = NO, 1 = YES # col11 = no. matings (only matters for males)# col12 = offspring produced: NA = NO, 1 = YES pop_matrix <-matrix(nrow = Starting_pop_size*2, # pop expands with initial repro pulsencol =12)# ID & Family ID pop_matrix[1:Starting_pop_size, 1:2] <-1:Starting_pop_size# assign sex pop_matrix[1:Starting_pop_size, 3] <-rbinom(n = Starting_pop_size, 1, prob =0.5)# female_starting_genotype pop_matrix[pop_matrix[,3] <1&!is.na(pop_matrix[,3]), 4] <-0# male_starting_genotype pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,3]), 4] <-0# assign mortality rates pop_matrix[1:Starting_pop_size, 5] <-1/baseline_mean_lifespan# set the unused rows to state -Inf pop_matrix[(Starting_pop_size +1):nrow(pop_matrix), 7] <--Inf# mate count pop_matrix[1:Starting_pop_size, 11] <-0# offspring production status pop_matrix[1:Starting_pop_size, 12] <-NA# populate breeding sites# the starting no. of females generally exceeds the number of breeding sites, which is starting_pop_size/f. The code below selects the initial breeding site holdersif(nrow(pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,3]),]) > breeding_sites){ initial_breeders <-head(pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,3]),1], breeding_sites) pop_matrix[initial_breeders,10] <-1# take advantage of ID = row number for initial pop } else{pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,3]), 10] <-1}# Initialise counter for the results table next_update <-0# keep track of when to update the results next_row <-0# keep track of which row to update# Initialise the Individual_ID and Family_ID counters Individual_ID_counter <- Starting_pop_size Family_ID_counter <- Starting_pop_size # each individual descends from a distinct family at onset# Initialise the timer t t <-0# Set initial pop size and freq of I allele for results table Prop_I <-0 pop_size <- Starting_pop_size total_female_deaths <-0 mated_female_deaths <-0# Start population without the I allele to generate family structure# Flips to 1 at mutant intro time point mutant_introduced <-0 keep_going <-TRUE# if the inbreeding allele fixes or goes extinct, this will change to false and the while loop will quit early# With the initial population ready to go, start the timer and let the simulation run.while(t <= time_end & keep_going){#print(paste0("Population size = ", pop_size, # ", breeders = ", sum(pop_matrix[,10], na.rm = T), # ", time = ", round(t, 3), ", Prop I =", Prop_I, ", mutation events =", mutant_introduced))# find next event # next death: this is the sum of the mortality rates for all individuals in the population next_death <- t +rexp(n =1, rate =sum(pop_matrix[, 5], na.rm = T))# next receptive mating encounter# find no. of females in mating pool & separate by encounter experience receptive_females_first_encounter <- pop_matrix[pop_matrix[,3] >0&is.na(pop_matrix[,6]) &is.na(pop_matrix[,7]),, drop =FALSE] receptive_females_second_encounter <- pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,6]) &is.na(pop_matrix[,7]),, drop =FALSE]# find no. of males in mating pool receptive_males <- pop_matrix[pop_matrix[,3] <1&is.na(pop_matrix[,7]),, drop =FALSE]# Find the time the next encounter occurs: plug the sum of the rates into the exponential function. # The population level encounter rate is the product of the rate at which a single male finds a single female, the number of receptive females in the population, and the number of receptive males in the population next_first_encounter <- t +rexp(n =1, rate = v*nrow(receptive_females_first_encounter)*nrow(receptive_males)) next_secondary_encounter <- t +rexp(n =1, rate = v*nrow(receptive_females_second_encounter)*nrow(receptive_males))# time in - Inf, Inf and Na are possible options that the code can handle next_time_in <-min(pop_matrix[is.finite(pop_matrix[,7]),7])# find which event happens next and update t t <-pmin(next_death, next_time_in, next_first_encounter, next_secondary_encounter, next_update, # update the populationna.rm =TRUE) # ... if a rate is 0, NaN produced.if(t == next_update &!is.na(next_update)){# record time, I prop and pop size results_matrix[next_row+1,1] <- t results_matrix[next_row+1,2] <-round(Prop_I, 4) results_matrix[next_row+1,3] <- pop_size # popsize results_matrix[next_row+1,4] <-round(mated_female_deaths / total_female_deaths, 3) next_update <- next_update +0.25 next_row <- next_row +1 total_female_deaths <-0# reset the count mated_female_deaths <-0# reset the count }if(t == next_death){# remove an individual from the pop who_died <-sample_vec(size =1, # choose onex = pop_matrix[!is.na(pop_matrix[,1]),1], # subset to current popprob = pop_matrix[!is.na(pop_matrix[,5]),5]) # weight by mortality rate# add a death if it was a femaleif(nrow(pop_matrix[pop_matrix[,1] == who_died &!is.na(pop_matrix[,1]) & pop_matrix[,3] >0,, drop =FALSE]) >0){total_female_deaths <- total_female_deaths +1}# add virgin female deathsif(nrow(pop_matrix[pop_matrix[,1] == who_died &!is.na(pop_matrix[,1]) & pop_matrix[,3] >0&is.infinite(pop_matrix[,7]),, drop =FALSE]) >0){mated_female_deaths <- mated_female_deaths +1}# remove individual from pop matrix pop_matrix[pop_matrix[,1] == who_died, 7] <--Inf# NA means time-in here, so special edit required pop_matrix[pop_matrix[,1] == who_died, c(1:6, 8:12)] <-NA# re-order to make steps like adding offspring easier later on pop_matrix <- pop_matrix[order(pop_matrix[,1]),] }# check if there are free breeding sites and whether females are available to fill them current_breeders <-sum(pop_matrix[, 10], na.rm = T)# get list of IDs for floating females floating_females <- pop_matrix[!is.na(pop_matrix[,1]) &# alive pop_matrix[,3] >0&# femaleis.na(pop_matrix[,10]), # non-breeding1] # return the IDs only# If so, recruit a new breeder# All prospective females have equal probabilityif(current_breeders < breeding_sites &length(floating_females) >0){# assign the new breeder new_breeder <-sample_vec(size =1, # choose onex = floating_females) # subset to floaters pop_matrix[pop_matrix[,1] == new_breeder, 10] <-1 }if(t == next_time_in &!is.na(next_time_in)){ # a male re-enters the mating pool pop_matrix[pop_matrix[,7] == next_time_in, 7] <-NA# change to receptive }#### matingif(t == next_first_encounter &!is.na(next_first_encounter)){# does first encounter lead to (inbred) mating?# Determine whether a heterozygote inbreeds on this occasion. # Depends on genotype if this matters heterozygote_inbreeds <-rbinom(1, 1, prob = dominance)# which female female_ID <-sample_vec(receptive_females_first_encounter[,1], 1)# get meta-data female <-subset(pop_matrix, pop_matrix[,1] == female_ID)# how many inbreeding alleles does she carry? alleles_female <- female[,4] mates <-NULL# reset this every time as a safeguard - MAYBE REMOVE?# find brothers that are in the mating pool brothers <- pop_matrix[pop_matrix[,2] == female[, 2] &# find family members pop_matrix[,3] <1&# that are maleis.na(pop_matrix[,7]) &# and in the mating pool!is.na(pop_matrix[,1]), # remove NAs1] # find the specific brother - if there aren't any, inbreeding does not happenif(length(brothers) >0){# choose brother randomly chosen_brother <-subset(pop_matrix, pop_matrix[,1] ==sample_vec(size =1, x = brothers))# how many inbreeding alleles does he carry? alleles_brother <- chosen_brother[,4] brother_ID <- chosen_brother[,1] }else{alleles_brother <-0} # we need this for the next if statement# now determine whether inbreeding occurs:# which individual expresses the allele# does that individual have the allele# is it expressed (depends on genomic region, no. copies and dominance)if(# female expression determines outcome# dominance doesn't matterlength(brothers) >0& sex_expressed >0& female_max_I == alleles_female |# dominance matterslength(brothers) >0& sex_expressed >0&0< alleles_female & alleles_female < female_max_I & heterozygote_inbreeds >0|# male expression determines outcome# dominance doesn't matterlength(brothers) >0& sex_expressed <1& male_max_I == alleles_brother |# dominance matterslength(brothers) >0& sex_expressed <1&0< alleles_brother & alleles_brother < male_max_I & heterozygote_inbreeds >0){# do inbreeding# update the pop matrix# female pop_matrix[pop_matrix[,1] == female_ID, 6] <-1# relative has been encountered pop_matrix[pop_matrix[,1] == female_ID, 7] <-Inf# female leaves mating pool pop_matrix[pop_matrix[,1] == female_ID, 8] <-1# inbreeding occurs pop_matrix[pop_matrix[,1] == female_ID, 9] <- alleles_brother # mates genotype# male pop_matrix[pop_matrix[,1] == brother_ID, 7] <- t + refractory_period # male leaves mating pool pop_matrix[pop_matrix[,1] == brother_ID, 8] <-1# inbreeding occurs pop_matrix[pop_matrix[,1] == brother_ID &!is.na(pop_matrix[,1]), 11] <- pop_matrix[pop_matrix[,1] == brother_ID &!is.na(pop_matrix[,1]), 11] +1 } else{# inbreeding is avoided# females that had no receptive brother to encounter are recorded as having had their chance for inbreeding early in life. When the male refractory period != 0, this is possible but unlikely (because all siblings are produced at the same time). Most commonly, this will occur when a female produces an all-female brood (0.03125 probability when f=5) pop_matrix[pop_matrix[,1] == female_ID, 6] <-1# relative has been encountered } }if(t == next_secondary_encounter &!is.na(next_secondary_encounter)){ # If the individual has already encountered a sibling, don't swap and let encounter proceed. # which female female_ID <-sample_vec(receptive_females_second_encounter[,1], 1)# get meta-data female <-subset(pop_matrix, pop_matrix[,1] == female_ID)# how many inbreeding alleles does she carry? alleles_female <- female[,4] # which male male_ID <-sample_vec(receptive_males[,1], 1)# get meta-data male <-subset(pop_matrix, pop_matrix[,1] == male_ID)# how many inbreeding alleles does he carry? alleles_male <- male[,4] # If the pair happen to be siblings, check if they inbreed # Determine whether a heterozygote inbreeds on this occasion. # Depends on genotype if this matters heterozygote_inbreeds <-rbinom(1, 1, prob = dominance)if(# female expression determines outcome# dominance doesn't matter female[,2] == male[,2] & sex_expressed >0& female_max_I == alleles_female |# dominance matters female[,2] == male[,2] & sex_expressed >0&0< alleles_female & alleles_female < female_max_I & heterozygote_inbreeds >0|# male expression determines outcome# dominance doesn't matter female[,2] == male[,2] & sex_expressed <1& male_max_I == alleles_male |# dominance matters female[,2] == male[,2] & sex_expressed <1&0< alleles_male & alleles_male < male_max_I & heterozygote_inbreeds >0){# do inbreeding# update the pop matrix# female pop_matrix[pop_matrix[,1] == female_ID, 7] <-Inf# female leaves mating pool pop_matrix[pop_matrix[,1] == female_ID, 8] <-1# inbreeding occurs pop_matrix[pop_matrix[,1] == female_ID, 9] <- alleles_male # mates genotype# male pop_matrix[pop_matrix[,1] == male_ID, 7] <- t + refractory_period # male leaves mating pool pop_matrix[pop_matrix[,1] == male_ID &!is.na(pop_matrix[,1]), 11] <- pop_matrix[pop_matrix[,1] == male_ID &!is.na(pop_matrix[,1]), 11] +1 } else{# do outbreeding# update the pop matrix# female pop_matrix[pop_matrix[,1] == female_ID, 7] <-Inf# female leaves mating pool pop_matrix[pop_matrix[,1] == female_ID, 9] <- alleles_male # mates genotype# male pop_matrix[pop_matrix[,1] == male_ID, 7] <- t + refractory_period # male leaves mating pool pop_matrix[pop_matrix[,1] == male_ID &!is.na(pop_matrix[,1]), 11] <- pop_matrix[pop_matrix[,1] == male_ID &!is.na(pop_matrix[,1]), 11] +1 } }# Consequences of death and mating: reproduction# check if a female can now produce offspring, either because they're previously mated and have secured a breeding site or because they already hold a breeding site and have now mated# make sure that previous breeders are excluded new_mated_breeder <- pop_matrix[is.infinite(pop_matrix[,7]) &# mated!is.na(pop_matrix[,10]) &# holds breeding siteis.na(pop_matrix[,12]),, drop =FALSE] # hasn't reproducedif(nrow(new_mated_breeder) >0){# add offspring to the population# each mated female that holds a breeding site produces f offspring# first check whether the mutant I allele should be addedif(mutant_introduced < mutation_events & t > mutation_time){ which_sex <-rbinom(1, 1, prob =0.5)if(chromosome =="A"& which_sex ==1| chromosome =="X"& which_sex ==1| chromosome =="Z"& which_sex ==1){ new_mated_breeder[4] <-1 }if(chromosome =="A"& which_sex ==0| chromosome =="X"& which_sex ==0| chromosome =="Z"& which_sex ==0){ new_mated_breeder[9] <-1 }if(chromosome =="W"| chromosome =="C"){ new_mated_breeder[4] <-1 }if(chromosome =="Y"| chromosome =="P"){ new_mated_breeder[9] <-1 } mutant_introduced <- mutant_introduced +1 } next_row_to_fill <-length(pop_matrix[!is.na(pop_matrix[,1]),1]) +1 last_row_to_fill <- next_row_to_fill + f -1 next_ID <- Individual_ID_counter +1 last_ID <- Individual_ID_counter + f Family_ID_counter <- Family_ID_counter +1# assign IDs pop_matrix[next_row_to_fill:last_row_to_fill, 1] <- next_ID:last_ID# assign all offspring to a single family pop_matrix[next_row_to_fill:last_row_to_fill, 2] <- Family_ID_counter# assign sex and genotype using our mating table sampling function offspring_genos <-sample_mating_table(inheritance_scheme, f, mother = new_mated_breeder) pop_matrix[next_row_to_fill:last_row_to_fill, 3] <- offspring_genos[,2] pop_matrix[next_row_to_fill:last_row_to_fill, 4] <- offspring_genos[,1]# assign mortality ratesif(is.na(new_mated_breeder[8])){ pop_matrix[next_row_to_fill:last_row_to_fill, 5] <-1/baseline_mean_lifespan } else{ # apply effect of inbreeding depression pop_matrix[next_row_to_fill:last_row_to_fill, 5] <-1/(baseline_mean_lifespan + D) }# fill in the mating and breeding site details - everyone starts as a floating virgin pop_matrix[next_row_to_fill:last_row_to_fill, 6:10] <-NA# mate count pop_matrix[next_row_to_fill:last_row_to_fill, 11] <-0# update the mothers offspring production status pop_matrix[pop_matrix[,1] == new_mated_breeder[1], 12] <-1# update the individual ID counter (redundant but more readable to do this here) Individual_ID_counter <- last_ID } # Calculate the frequency of the I allele, quit early if I fixes or goes extinct pop_size <-nrow(pop_matrix[!is.na(pop_matrix[,1]),, drop =FALSE]) # use this to update the results n_females <-nrow(pop_matrix[!is.na(pop_matrix[,1]) & pop_matrix[,3] >0,, drop =FALSE]) n_males <- pop_size - n_females# calc allele freq if autosomal locus if(chromosome =="A"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/(pop_size*2) # x2 because diploid }# calc allele freq if W locus if(chromosome =="W"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/n_females }# calc allele freq if Y locus if(chromosome =="Y"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/n_males }# calc allele freq if X locus if(chromosome =="X"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/(n_females*2+ n_males) }# calc allele freq if Z locus if(chromosome =="Z"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/(n_females + n_males*2) }# calc allele freq if C locus if(chromosome =="C"| chromosome =="P"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/pop_size }# quit conditionif(mutant_introduced >0& Prop_I >0.9| mutant_introduced >0& Prop_I ==0| pop_size <2){keep_going <-FALSE} } results_matrix[next_row+1,1] <- t results_matrix[next_row+1,2] <-round(Prop_I, 4) results_matrix[next_row+1,3] <- pop_size results_matrix[next_row+1,4] <-round(mated_female_deaths / total_female_deaths, 3) results_matrix <- results_matrix[-(next_row+2:nrow(results_matrix)),]# save results as a csv. results_matrixwrite.csv(results_matrix,paste("results/rowID_", parameter_space_ID, chromosome, ".csv", sep =""))#write.csv(results_matrix,# paste("sim_results/rowID_", # parameter_space_ID, # chromosome, ".csv", # sep = ""))}
Run the simulation
In practice, we ran the simulations on JGU’s Mogon computing cluster. See the HPC_inbreeding_script.R and the batch script run_inbreeding_sim. To run the simulation for a single parameter space, you could run continuous_time_simulation(1, parameters_P, offspring_genotypes_P)
Load the results
Build a helper function for extracting results
Code
# this function loads the individual runs and joins them into a single tibblefiles <-list.files(path ="sim_results") %>%str_subset("P") # change this to load the desired filesresults_reader <-function(x){read_csv(paste0("sim_results/", x)) %>%mutate(parameter_space_ID = x)}
First, let’s load the results for the coarse parameter space.
The number of individuals present in the population for each simulation run is plotted below. We include this primarily as a sanity check. Browse the tabs to view population trajectories when inbreeding acceptance alleles occur at different genomic locations. The dotted line indicates the appearance of the inbreeding acceptance alleles in the population. Each run ends when either i) the inbreeding acceptance allele is purged, ii) the inbreeding acceptance allele approaches fixation or iii) the population goes extinct.
Figure S1. The effect of mate searching efficiency on female matelessness. Each line shows the outcome of a single simulation run across the first 50 time units. Runs differ in parameter values. Orange and blue trajectories indicate simulation runs where males have low and high mate searching capabilities, respectively. The dotted line indicates the appearance of the inbreeding acceptance alleles in the population.
Building Figure 2
Prep data for plotting. We first extract the final \(I\) allele frequencies from each simulation run and categorise the outcome as either i) Population extinction, ii) purging of the \(I\) allele or iii) \(I\) approaches fixation.
Figure 2. Results generated by stochastic population genetic simulations. Panels show the conditions where dominant inbreeding alleles invade to a high frequency (0.9) when females (left column) or males (right column) control the outcome of mating interactions, for genes found on W, maternally inherited cytoplasmic, X, Autosomal, Z, paternally inherited cytoplasmic and Y chromosomes. Inbreeding depression is expressed as the mean proportion of lifespan lost relative to an outbred individual, while the male refractory period is measured as the proportion of an outbred individuals maximum lifespan that is spent out of the mating pool after mating. The thick dashed lines show \(\delta_{\mathrm{f}, i}^*\) and \(\delta_{\mathrm{m}, i}^*\) when inbreeding is sporadic, while the thin dashed line shows these thresholds at the frequent inbreeding extreme (if they differ). Note that each pixel reports the outcome of precisely two simulation runs, one where male-female encounters are common and another where they are rare. The mosaic nature of the plot reveals the stochasticity of outcomes when a small number of mutants tried to invade per run. While inbreeding is therefore never guaranteed to invade, the plot allows identifying parameter regions where it never spreads, as opposed to often doing so.
Figure S4. Results generated by stochastic population genetic simulations. Panels show the conditions where inbreeding invades to a high frequency (0.9) at autosomal loci, when females (left column) or males (right column) control the outcome of mating interactions, when inbreeding acceptance alleles are either completely recessive, additive or dominant. Inbreeding depression is expressed as the mean proportion of lifespan lost relative to an outbred individual, while the male refractory period is measured as the proportion of an outbred individuals maximum lifespan that is spent out of the mating pool after mating. The thick dashed lines show \(\delta_{\mathrm{f}, i}^*\) and \(\delta_{\mathrm{m}, i}^*\). Note that each pixel reports the outcome of precisely two simulation runs, one where male-female encounters are common and another where they are rare. As alleles become more recessive, they are less likely to be expressed while rare, making it more difficult to overcome the drift barrier.
Figure S5. Results generated by stochastic population genetic simulations. Panels show the conditions where inbreeding invades to a high frequency (0.9) at X-linked loci, when females (left column) or males (right column) control the outcome of mating interactions, when inbreeding acceptance alleles are either completely recessive, additive or dominant. White fill indicates dominance scenarios we do not explore, due to hemizygosity. All other details are as in Figure S2.
Figure S6. Results generated by stochastic population genetic simulations. Panels show the conditions where inbreeding invades to a high frequency (0.9) at Z-linked loci, when females (left column) or males (right column) control the outcome of mating interactions, when inbreeding acceptance alleles are either completely recessive, additive or dominant. White fill indicates dominance scenarios we do not explore, due to hemizygosity. All other details are as in Figure S2.
Plotting results over the fine-scale parameter space
Figure S2. Results generated by stochastic population genetic simulations, for levels of inbreeding depression where \(\delta_{\mathrm{f}, X}^*\) and \(\delta_{\mathrm{m}, X}^*\) are predicted to depend on the emergent values of \(F_\mathrm{f}\) and \(F_\mathrm{m}\). Panels show the conditions where inbreeding invades to a high frequency (0.9) – or in very rare cases should eventually reach this cut-off, but had not after 1000 time units – when females (left column) or males (right column) control the outcome of mating interactions, for genes found on X and autosomal chromosomes. The autosomal situations act as controls, with inbreeding acceptance predicted to invade across a larger fraction of the plotted space compared with the X-linked case, if frequent inbreeding affects the invasion dynamics. Note that these panels are zoomed in versions of the relevant panels shown in Figure 2, simulated with a finer resolution. All other details are as in Figure 2.
Figure S3. Results generated by stochastic population genetic simulations, for levels of inbreeding depression where \(\delta_{\mathrm{f}, Z}^*\) and \(\delta_{\mathrm{m}, Z}^*\) are predicted to depend on the emergent values of \(F_\mathrm{f}\) and \(F_\mathrm{m}\). Panels show the conditions where inbreeding invades to a high frequency (0.9) – or in cases that should eventually reach this cut-off, but had not after 1000 time units – when females (left column) or males (right column) control the outcome of mating interactions, for genes found on Z and autosomal chromosomes. The autosomal situations act as controls, with inbreeding acceptance predicted to invade across a smaller fraction of the plotted space compared with the Z-linked case, if frequent inbreeding affects the invasion dynamics. Note that these panels are zoomed in versions of the relevant panels shown in Figure 2, simulated with a finer resolution. All other details are as in Figure 2.
